Non-equilibrium transport response from equilibrium transport theory 
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We propose a simple scheme that describes accurately essential non-equilibrium effects in nanoscale electron- 
ics devices using equilibrium transport theory. The scheme, which is based on the alignment and dealignment 
of the junction molecular orbitals with the shifted Fermi levels of the electrodes, simplifies drastically the cal- 
culation of current-voltage characteristics compared to typical non-equilibrium algorithms. We probe that the 
scheme captures a number of non-trivial transport phenomena such as the negative differential resistance and 
rectification effects. It applies to those atomic-scale junctions whose relevant states for transport are spatially 
placed on the contact atoms or near the electrodes. 
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Nano-scale devices tend to suffer from a serious lack of re- 
producibility from one research group to another, and from 
a lack of durability and robustness [ ]J |2fl. However, the con- 
sistent and thorough work of a range of experimental groups 
have enabled to reach some consensus on the conductance 
values and variability of a few specific junctionsj3l-|5t]. In 
addition, the development of simulation codes based on a 
combination of density functional theory (DFT) and the non- 
equilibrium Green's function formalism fNEGF1ll5l- [l2ll have 
enabled to assist the above experiments with theoretical in- 
sights and predictions. However, the size and complexity 
of the experimentally active part of the junction is typically 
too large, and the above codes need to assume a number of 
simplification regarding the size, geometry, number of feasi- 
ble atomic arrangements and the electronic correlations. Re- 
cently, some of the above codes have been combined with 
classical force field programs, enabling the simulation of 
much larger systems and the sampling of a much wider set of 
atomic arrangements ITil [Till . Nevertheless, present current- 
voltage I — V characteristics are still computed using NEGF 
techniques, which are rather cumbersome and extremely com- 
puter hungry. In other words, the NEGF is a serious bottle- 
neck hindering the deployment of realistic transport simula- 
tions at the nanoscale. In contrast, equilibrium transport tech- 
niques are much simpler, better founded on physical grounds 
and computationally drastically less demanding [6|. How- 
ever, they are completely unable to reproduce current- voltage 
curves of nano-scale device s lfl5ll . This is so because these 
techniques assume that the electronic states and the Hamilto- 
nian at the junction do not change under the application of a 
voltage. However, a voltage bias drives a flow of electrons into 
and out of the junction, driving the device out of equilibrium 
and possibly even changing its quantum state and sometimes 
even its atomic arrangement. For example, the electric field 
originated by the voltage bias can also shift the energy posi- 
tion of the molecular levels by the Stark effect, depending on 
the polarizability of each molecular level or on whether the 
molecule has itself a polar nature 1 3 171. 



transport techniques that accounts for a large set of non- 
equilibrium transport phenomena such as negative differen- 
tial resistance (NDR)[ 18] or rectification effects. The perfor- 
mance of our method depends on the energy position, elec- 
trode coupling and spatial localization of the molecular or- 
bitals that are responsible for the transport properties. This 
works specially well for systems where the closest states to 
the Fermi level (HOMO or LUMO) are located at the contact 
atoms or near the electrodes. We discuss below the physi- 
cal mechanisms at work for several molecular-scale functions. 
Our scheme enables the simulation of I — V curves which 
compare favorably to those obtained using NEGF techniques. 
As a consequence, we also propose that a number of transport 
phenomena regarded up to now as genuine manifestations of 
non-equilibrium behavior of a junction could be recognized 
as quantitative but not qualitative modifications of its equilib- 
rium state. 

We use Caroli's scheme [20] and split the device into left 
and right electrodes, and a extended molecule. The extended 
molecule includes several principal layers of the electrodes, 
as well as the molecule itself. The full Hamiltonian H can be 
written as 
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We show here a suitable modification of equilibrium 
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where the matrices Hi, Hr, Him, Hrm and their corre- 
sponding overlap matrix blocks, indicate respectively the left- 
and right-hand side leads Hamiltonians and the coupling ma- 
trices between the leads and the extended molecule. The ma- 
trix elements of the extended molecule Hamiltonian H^ run 
over all the orbitals residing in the extended molecule region. 
They depend on the non-equilibrium density matrix p, which 
must be calculated self-consistently for each voltage V us- 
ing cumbersome and computationally very demanding NEGF 
techniquesS&tl2l. 



Typeset by REVTbX 



2 



< 

a 

1) 

H 
3 

u 



=L 

fcs 
u 



- — Non Equilibrium 
■ ■ ■ Equilibrium 
On-site shift 




1 


1 1 1 ' 

- (a2) 


(al) 












1 1 h 

(bl) 

* '^^^ 




i(b2) J 


',1,1 


1 





1 2 3 1 
Voltage (V) 



2 3 



90 

60 

30 


8 
6 
4 
2 




FIG. 1 : (Color online) / — V curves of a junction made of two car- 
bon chains, each terminated by a nitrogen (a) or a sulfur atom (b). 
These two atoms are separated by a physical gap of 3.7 A. Panels 
(al) and (bl) correspond to fully relaxed coordinates. Panels (a2) 
and (b2) correspond to chains where the terminating nitrogen and 
sulfur atoms have been separated by 0.5 A away from their relaxed 
positions. Solid and dotted lines correspond to NEGF and equilib- 
rium simulations, respectively. Dashed lines correspond to simula- 
tions using the scheme proposed in this article. 



We propose here that in many situations p and its de- 
rived Hm[p] need only be computed at zero voltage, and 
that the effect of a finite bias can be accounted for by the 
alignment or dealignment of the molecular orbitals at the ex- 
tended molecule region with the energy levels of the elec- 
trodes (which are shifted by ±eV/2 by the bias voltage). 
Mathematically, we apply a simple shift to the Hamiltonian 
matrix elements from Hm to "Hm + eV% Sm> where the lo- 
cal shifts Vi depend on the junction nature. For example, in 
a highly-transparent junction, the shifts can be modeled by 
a linear voltage ramp connecting the matrix elements of the 
orbitals at the outest layers. In contrast, most nanoscale junc- 
tions contain a molecule that is connected to the left and right 
electrodes by weaker links, which means the details of the 
coupling of the the electrodes to the molecular orbitals re- 
sponsible for the transport must be taken into account. Our 
scheme is specially suited for junctions where the HOMO and 
LUMO orbitals are spatially localized at the linker atoms join- 
ing the molecule with the electrodes. In this case, we take 
Vi = ±e V/2 for all orbitals starting at the outest layers in 
the extended molecule and up to the linker atoms, and Vi — 
for all the remaining inner orbitals. We show now how these 
shifts allow to mimic accurately a number of simple and not 
so simple junctions displaying non-trivial negative differential 
resistance, as well as rectification effects for some asymmetric 
molecules. 

Notice that in NEGF simulations the left and right elec- 
trodes are defined as semi-infinite electron reservoirs so that a 
steady electron flow is established. The chemical potentials at 
the left and right electrodes are shifted by ±eV/2 and a per- 
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FIG. 2: (Color online) Transmission curves at different voltages for 
the carbon junctions shown in panels (a2) and (b2) of Fig. (T}. Panels 
(al) and (bl) correspond to a full NEGF calculation. Panels (a2) and 
(b2) correspond to the scheme proposed in this article. 



manent potential drop is produced somewhere in the extended 
molecule. These features are reproduced by our scheme. We 
have also checked that they are not reproduced by applying 
an electric field along the transport direction in equilibrium 
simulations of finite-sized systems that represent the extended 
molecule. This is so because in these later cases the junction 
has a finite length and equilibrium conditions are reestablished 
at the end of the simulation. We have verified that the ap- 
plied electric field shifts initially the electronic states as in our 
scheme, but because the whole device has a common Fermi 
energy, a positive charge surplus accumulates at the negative 
electrode and vice versa, thereby screening the potential cre- 
ated by the electric field and moving back the energy levels 
towards their original positions. 

We discuss first two simple examples to illustrate the 
scheme. These systems are made of two semi-infinite carbon 
chains terminated each by either a sulfur or a nitrogen atom 
and separated initially by a vacuum region of 3.7 A. We apply 
a shift of ±e V/2 to each Hamiltonian matrix element at the 
extended molecule region up to and including the terminating 
atoms as proposed above. The left panels in Fig. (Q~|) show 
the I — V characteristics of the two junctions, which display a 
conventional ohmic behavior. However, these junctions show 
NDR when the terminating atoms are slightly displaced away 
from the chains[21]. The I — V curves of the two junctions 
are plotted in the right panels of Fig. (Q]). The NDR response 
shown is originated by a non-trivial energy shift and splitting 
of the resonances appearing in the transmission coefficients 
T(E) as the voltage bias is ramped up. This non-trivial evolu- 
tion of the transmission resonances is shown in the left panels 
of Fig. (0 and could be thought of a genuine non-equilibrium 
effect at first sight. That this is so could also be inferred from 
the / — V curves computed from the equilibrium equation 
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where T(E) is calculated at zero voltage. Such I — V curves 
are indeed quite different from those computed using stan- 
dard NEGF techniques, as shown in the right panels of Fig. 
(0. However, the figure also shows that Eq. (0 can produce 
extremely accurate results if T(E) is calculated using a mod- 
ified Hamiltonian Hm that follows our scheme. Furthermore, 
the right panels in Fig. (0 show that the non-trivial bias evo- 
lution of the transmission resonances is also followed by our 
scheme with high fidelity. As a consequence, we propose that 
the resonance splitting and shift shown in the figure is only 
apparent. What actually happens is that the single resonance 
is cut-off by an opposite shift of the on-site energy levels at 
opposite sides of the junctions, and is therefore not a genuine 
non-equilibrium effect. 

The scheme proposed cannot be applied to all physical sit- 
uations. An example where the scheme fails occurs when 
the molecular orbital responsible for transport (HOMO or 
LUMO) is well localized in the central region of the molecule. 
In this case, the resonance associated with that molecular or- 
bital could be pinned to the Fermi level of either the left or 
the right electrode 112211 - However, even in this scenario the 
scheme proposed improves slightly the /—^curves compared 
to raw equilibrium simulations, because it takes into account 
a non-equilibrium mechanism that is present in all nano-scale 
junctions. We discuss next a series molecular junctions where 
different polyyne derivative are contacted by (001) gold elec- 
trodes to illustrate this argument. Results for the synthesis of 
the molecules and the transport properties of the devices have 
been reported in Ref . B23I1 . 

Fig. (0 (a) shows the I — V curves of a pyridine-polyyne- 
pyridine junction, where the polyyne has n = 3 carbon pairs. 
Following our scheme, we apply a rigid shift to the orbitals 
of the extended molecule up to and including the nitrogen 
atoms at the two pyridine rings. The resulting I — V curve, 
plotted in Fig. (0 (a), shows that the scheme fails calami- 
tously to describe the NDR response of this junction, even 
though it improves the equilibrium-based calculation of the 
I — V curve, which always increases. Furthermore, we have 
tested that extending the shift up to the inner carbon atoms 
in the pyridine rings does not improve matters. In contrast, 
the scheme provides reasonably accurate curves for similar 
polyyne-based junctions where the pyridines are replaced by 
sulfur atoms (Fig. (0 (b)) or by sulfur-benzene rings (Figs. 
(c-d)). The corresponding T(E) curves for the pyridine- 
polyyne-pyridine and sulfur-benzene-polyyne-benzene-sulfur 
junctions, calculated with NEGF and with our scheme, are 
shown in Fig. (0. This figure shows how in the first case, a 
sharp resonance exists, which moves to higher energies as the 
bias increases while maintaining its shape. In contrast, our 
scheme predicts that the resonance does not shift with voltage 
because it is placed at the center of the molecule, as we show 
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FIG. 3: (Color online) I — V curves of gold (001) junctions 
sandwiched by pyridine-polyyne-pyridine (a), sulfur-polyyne-sulfur 
(b), sulfur-benzene-polyyne-sulfur (c) and sulfur-benzene-polyyne- 
benzene-sulfur (d) molecular wires, where the polyyne has n — 3 
carbon pairs. The molecules are sketched above or below each cor- 
responding graphic, together with density contours that highlight the 
localization of the relevant molecular orbital. Solid lines correspond 
to non-equilibrium calculations, while dashed lines correspond to 
calculations performed with the scheme proposed in this article. The 
on-site energy levels have been shifted up to the atoms indicated by 
the arrows. 



in Fig. (0 (a2). The impact of a finite bias in the second junc- 
tion, however, is manifested in the destruction of the HOMO 
resonance, shown in Fig. (0 (b), which is now well captured 
by our scheme. 

The criterion to discern when the proposed scheme can be 
used is related to the nature and features of the molecular level 
responsible for the transport properties. Those are its spatial 
placement in the molecule, its energy positioning as well as its 
coupling to the electrodes. For the pyridine-polyyne-pyridine 
molecule, the LUMO orbital is made essentially of carbon it 
states and is located on the central part of the molecule. When 
the bias is applied, the voltage drop does not occur at this or- 
bital because the highly conjugated character of the molecule 
screens the potential. Therefore the applied bias has a small 
impact in the LUMO orbital and as a result the shape of the 
associated transmission resonance does not change. The res- 
onance shifts up in energy because there is a net transfer of 
charge from the electrodes towards the molecule. This moves 
all levels upwards and pins the LUMO orbital to the nega- 
tive electrode. For the sulfur-benzene-polyyne-benzene-sulfur 
junction, the relevant orbital is the HOMO, which is mostly 
located at the sulfur atoms. The atomic sulfur states couple 
across the molecule and produce bonding and anti-bonding 
states. When the bias is ramped up, those sulfur states fol- 
low the chemical potential of each electrode because they are 



/I 



strongly coupled to the leads. As a consequence, the bonding 
and anti-bonding molecular levels are destroyed. This is seen 
in the transmission as a reduction of the HOMO-related res- 
onance, which produces a NDR response. We note that this 
NDR behavior does not occur if the molecule does not con- 
tain the benzene rings as shown in Fig. (0 (b). In this case 
the coupling between the sulfur states and the carbon chain 
is much stronger. This increases the hybridization between 
them, keeps the integrity of the bonding and anti-bonding 
states and decreases the impact of the shifts. An intermediate 
case between these two junctions corresponds to the asymmet- 
ric sulfur-benzene-polyyne-sulfur molecule shown in Fig. ® 
(c). This junction gives different I— V characteristics for each 
bias polarity, producing a rectifying response. In this case the 
relevant orbital is again the HOMO state, which is mostly lo- 
calized in the right sulfur atom and the carbon chain. The 
chain states are left unaffected when we ramp up the bias and 
apply the scheme to the sulfur on the right side and the ben- 
zene carbon atom closer to the chain. The scheme reproduces 
quite well both the NDR and the rectification effects. 

The scheme proposed is specifically suited for junctions 
where the states are well localized at the linking atoms or 
near the electrodes, because then the frontier orbitals are 
strongly affected by the applied bias. Large junctions are diffi- 
cult to simulate using NEGF techniques because they require 
huge amounts of computer resources and have frequently con- 
vergence problems. Our scheme simplifies the calculation 
of I — V curves drastically, because it relies on equilib- 
rium techniques, which are much simpler and computation- 
ally much less demanding. We have indeed verified that the 
I — V curves corresponding to a series of graphene junctions 
showing NDR response ll2lll are accurately described with our 
scheme. The recipe proposed in this article can be used when- 
ever the ground state at the molecule is not changed by the 
non-equilibrium conditions. In the electrostatic spin effect, for 
example, the spin multiplicity of the molecule is changed by 
the application of a finite bias[ 17], and as a result our scheme 
cannot be applied. 

In summary, we propose in this article that current-voltage 
characteristics of a wide array of electronic junctions can be 
accurately described by a scheme that combines conventional 
equilibrium transport techniques with suitable shifts of en- 
ergy levels up to the weak links existing in those junctions. 
This scheme can have drastic practical effects because equilib- 
rium codes are much simpler and computationally more effi- 
cient than non-equilibrium codes. We have corroborated very 
large computational speed-ups in the above examples. As a 
consequence, our recipe opens a route towards realistic trans- 
port simulations comprising much bigger or complex systems 
than those that could be feasibly simulated with current NEGF 
techniques. 
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FIG. 4: (Color online). Transmission curves at different voltages 
for the pyridine-polyyne-pyridine (top panels) and sulfur-benzene- 
poly yne-benzene- sulfur (bottom panels) junctions shown in Panels 
(a) and (c) of Fig. l[3j. Panels (al) and (bl) correspond to a full 
non-equilibrium calculation. Panels (a2) and (b2) correspond to the 
scheme proposed in this article. 
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